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ABSTRACT 


A program is developed to ascertain the response of a 
frictionless water column to internal waves. Using density 
strata aS input, the program selects ten frequencies at 
equal intervals within a spectrum of internal waves bounded 
above by the maximum Vaisala frequency and below by the 
inertial frequency. Ray paths within the column are plotted 
for these frequencies. The first ten normal modes for each 


frequency are computed. The first three modes are plotted. 
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Wee LNT RODUCTION 


Internal waves exist in all oceans, most bays and lakes; 
they vary widely in amplitude, period, and depth (La Fond, 
1962). Recently, advances have been made in interpreting the 
major features of the oceanic internal wave field (Wunsch and 
Mellen, 1970). 

Internal waves can be observed as density perturbations 
through a stable, stratified water column. Since temperature 
correlates with density in the majority of the upper ocean 
areas, internal waves often appear as oscillations about the 
region of the thermocline, exhibiting greatest amplitude and 
pisces pecan Mere No tueitmocline iS the most intens@. ‘the 
@scillation initiated within the water column is described in 
pemmoanoie the Vaisala or stability frecuency. This frequency 
may be derived as follows (Zalkan, 1966). 

Consider the force acting on a water parcel raised 


vertically from its equilibrium position: 


et gs 20° fa 
=¢ Q- (+57 Az) Jav-P 5 (Nae (1) 


where: g=gravitational acceleration, 
O = density, 
Pea veorurcal displacement , 
dV =minute volume of the parcel, and 
c=sound velocity. The first term on the right hand 
sadesis the restoring force due to buoyancy. The second term 


is tne contributing force of compressibility. The equation 


T 





for the induced vertical motion is: 
2 
a ZaVv “eS ZaV (2) 


The harmonic solution to equation (2) is the Vaisala frequency 


defined as 


(3) 





For oceanographic applications the water is effectively in- 


compressible and the second term under the radical is 


- (2% 
N(Z) = = Sy (4) 


The Vaisala frequency is greatest in magnitude where the 


negligible. Thus, 


thermocline is most intense (see Figure 1). 


itemwemrcal NOukoON iniwsated by the intern 
ing through the water column is governed by the differential 


equation (Fjeldstad, 1933): 
oan k- pe DD 
7(O <3) +> (N“ -§")W=0 (5) 
0 Z C a2 if = ee 


where: Wz=vertical water velocity, 

k=internal wave number, 

( = internal wave frequency, and 

f=inertial frequency. Equation (5) is a hyperbolic 
wave equation whose characteristics describe the direction 
of propagation of internal wave energy. With the appropriate 
boundary conditions that no vertical flow exists at the sur- 
face or at the bottom, the internal wave equation satisfies 


the criteria for a Sturmian system. The imposed boundary 
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conditions are satisfied by the simultaneous propagation of 
two waves having the same horizontal wave number, one with 
its vertical component of phase velocity up, the other with 
its vertical component down. The stipulation that these 
waves cancel at the boundaries leads to the structure of 
normal modes. These modes constitute a discrete set of wave 
numbers that are the eigenvalues of the Sturmian system. 

This thesis, using as a basis Fjeldstad's classic paper 
(Fjeldstad, 1933), developed a numerical program that assessed 
the response by a frictionless, incompressible water column 
to internal waves in a frequency band bounded above by the 
maximum Vaisala frequency and below by the inertial frequency. 
Station density strata was used as input. Frequencies were 
selected within the band and their ray paths were drawn. 
Also, the governing differential equation was numerically 
integrated and the resulting endpoints were refined to deter- 
mine the first ten normal modes. The first three modes were 
plotted. 

The program was first tested using density strata with 
known output. The program then treated data taken from ocean 
Stavulons o1f the coast of Monterey, California in mid-August, 


iPS EO. 
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ee Co acOhu 


A. DEVELOPMENT OF THE INTERNAL WAVE EQUATION 
For a frictionless, incompressible water column the 


governing equations (following Feldstad, 1933) are: 


Momentum 
OU ov 48 SE =o (6) 
oY 4 1u+5 SP -o (7) 
OW 10P 
gt FO Oz = 0 (8) 
Continuity 


— +3 +55 =0 (9) 


Density invariant 


Ku Biv hiv 8-0 (10) 


The X-axis is directed in the direction of wave travel, 
the Y-axis is in the horizontal transverse plane of travel, 
PacbettemZ-axte ts vertacal, positive upward. U, V, and W are 
the velocity components of the water in the X, Y, and Z axes, 
respectively. f is the coriolis parameter. FP is the pressure. 
Mre depth is assumed to be uniform. Z=0 is the bottom; Z=h 
is the undisturbed upper surface. For a steady state 


CoOndL Vion , 


dqU_ dV_ dW _ 
dima din (11) 
Einise. 
dP = -egdZ. (12) 


ile 





P can be represented as 
P=p ote f par, Gis 


where Po is the undisturbed pressure at the surface. Let 
© =0,(z) + © (x,¥,2,+), (14) 
where a is the undisturbed density and c, the perturbation. 


Therefore, pressure for the perturbation is: 


Pay ta) Gaze, (15) 


where 


= J, @a2 (16) 


Envoking a general assumption that velocities and perturbations 
of density and pressure are small enough that their products 


and squares may be neglected, equation (10) tay be written: 


i 5te,02'"° (17) 
tt onmiiromiibernal wave £orm, harmonic oscillations 
such as the following are considered: 
uj€, ,P, & cos @t-kx) exp (ay) 
V,WA sin (ft — kx) exp (ay) 
VietesLne "cCOnstany a controls the transverse envelope and 
allows for the possibility of edge waves. The sign of a 1s 


such that the wave amplitude decays away from the boundary. 


The governing equations now become: 


al 
(y U-fV+k—=0 (18) 
So 
4 2 (19) 
V+fU+ an-=0 19 
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W+tgt+5- y= 0 20 

‘ef Ro Q Zz ( ) 
dW 

kU+aVt+—7=0 (21) 


P 
Ken, ae 
U=— DD 
Wy 2 o8e Rp ae 
= 
k 
oO 5 (23) 


(24) 


where N =o" ra as previously defined. When one inserts 


es roan on 


2% Cae =) dZ, 
When P, is eliminated from (24) and (25) and the transverse 


boundaries are disregarded (i.e., set a=0), 


d-W aN. a5 i 


“3 (@- az az" of ally" -§*Jw=0 (26) 


Equation (26) is the internal wave equation describing the 
response by a water column to a plane internal wave of known 
frequency,0. The order of magnitude of the second term is 
small compared to the first and third. However, the numerical 
program can include the second term without any difficulty. 


ittiewmli be so retained. 





BeOS eELECATION OF RAY THEORY 
An analytical solution to equation (26) may be made by 


the transformation: 


Q 
W= (Z) exp Feb. (27) 


When the effect of the second term in eq. (26) is small, the 


analytical solution to the internal wave equation can, with 


a high degree of accuracy, be represented Bee, the 


Beewllavory portion. The oscillatory solution is: 
\° kK 2 
zz, * | D 5) (N -J \W=0 (28) 
g“-f 
Solutions to equation (28) are of the form: 
Y «Coe 
ee) 1 (kK, X+k, ye ee Y+q¢t) (29) 


where K, = horizonta1 wave number. 
Define © as the angle between the horizontal and the 


phase velocity. Then, 


k 2 SE ey) 2: 
008 @ =——4, 575 = S55 (30) 
2 2 No = Net 
Since 
ye Ly 2 k 
(2°) = (N* -9*) San (31) 
(x +k, ) 
% 2Ne 4 kA 2 
: Ky + ky* = 
where / 
os. ean eye: 
N 
fe, * + %,7)'/* (33) 





Figure 2 illustrates the angle between the horizontal 


and the phase velocity. The phase velocity, broken into X 


and Z components, is: 


ky Ky 


Uae, Os 
C RONEES = D1... 372 [kk N° +k, f | 


il 2 


and 


k 
2 D2 2 
Si = ey A [k, “N° +k.°f 
7. (1,2 +2) | 2 


Per the group velocity, 
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C, a5 
where 
2) Exc: ane 
Vic= (Sig i * Ok; J ) 
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~> AY, No ke : 
1 Zz 


The X-component of group velocity is: 
\_ ,2nt 44,24? 1/2 
ee (>) 
X~ gk) 2 eZ 


Ky +k, 
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Mage) No _ f° 





Kk) N +k, 
The Z-component of group velocity is: 
MeN ces fe 1/2 
2 hoe 16 


Define — as the angle made with the horizontal 


group velocity vector. Then, 


ties 
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je||2 Dee Ze Ny: 2 


(35) 


(39) 


(40) 


by the 








PIGURE 2. Illustration of Angle between 
Horizontal and Phase Velocity 
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tan ©, = (gy, --K, = Se (41) 


The ray trace part of the program drew the ray paths as 
depicted in equation (41) for ten frequencies selected within 
the inertio-gravity band. 


Determining the dot product of the phase and group 





moLocities: 
—> 
Cac) Cia 
g 
OE, 2 <2 2 
05 eerie f= ) k,k, n2 _ £2 , 
lic |? lx]? (k,°N* +k,°£*)? 2 
k, (ik, “NA a yg) 1/2 k, “k, (Ne _ 7°) 
= ) ( 3 Oe DES 173) 
lg] CS Naat) 


Thus the phase velocity is perpendicular to the group 
velocity. Also, 

k(ky° - ky) £° 
A Se aT ETSY LTE ESTE (43) 


The phase velocity and group velocity then have a relationship 


as seen in Figure 3. 


C. NORMAL MODE APPLICATION 

With the use of suitable boundary conditions, equation 
(26) may be solved to obtain modal structures. The applicable 
Conainlons are: 


Ve 0 ale ie —20) (44) 
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Figure 3. Relation between Phase and Group Velocities 
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and 


Ps 
dW k 


To a high degree of approximation the vertical velocity shear 
is negligible at the surface and the second condition becomes: 
W=-0 at Lah (45) 
Equations (26), (44), and (45) compose a Sturmian system. 
W, the vertical velocity, is the eigenfunction. For a given 
frequency, discrete wave numbers or eigenvalues are determined. 
The program utilized a fourth-order Runge-Kutta numerical 
method to obtain an end value after integration of the 
internal wave equation over the water column. The endpoint 


was refined using a Newton-Raphson iteration procedure. 


Ey, 





III. DEVELOPMENT OF THE MODEL 


A. MAIN PROGRAM 

The main program prints the station density stratifica- 
tion. The program then calls upon the three sub-programs to 
interpolate density between station depths, to plot ray paths 
for selected internal wave frequencies, and to print out and 


plot modal structures for these frequencies. 


B. SUBROUTINE INTERP 

pubroutime INTERP interpolates the density stratification 
for 400 incremental depths through the water column. The 
interpolation is accomplished by the use of a third-degree 
polynomial that matches points and slopes at each depth where 
data was obtained from a Nansen bottle. To guarantee density 
stability, the smaller of the slope above or below each point 
is chosenas the slope at that point. From the interpolation 
polynomial, INTERP computes the Vaisala frequency for each 


density value. 


C. SUBROUTINE RAY 

Subroutine RAY plots the interpolated density and the 
Vaisala frequency as a function of depth. RAY then determines 
the inertio-gravity frequencies at rates of 1/100, 1/10, 2/10, 
oy erm Oem >, Tro, 10, 7/10, 8/10, and 9/10 of the frequency 
band. The program then plots the characteristic paths over 


the depth of the column for each of these frequencies. 


20 





Additionally, subroutine RAY determines the percentage of 
the water column depth where the selected frequency exceeds 
the Vaisala frequency. This percentage represents that por- 
tion of the column where the Sturmian solution to the internal 
wave equation is exponential. <A greater percentage implies 
miereasing ditficulty in applying boundary conditions during 
the numerical integration. The percentage figure for each 


frequency is listed beneath the ray path plot. 


D. SUBROUTINE RUNK 

Subroutine RUNK utilizes a fourth-order Runge-Kutta 
numerical method to integrate the internal wave equation over 
the water column. The endpoint of the integration is refined 
by a Newton-Raphson iteration method until the endpoint is 
1/200 of the maximum value of the eigenfunction, W, determined 
during the integration scheme. For each selected frequency, 


piowterst approximation for the initial eigenvalue 1s 


v 


Se (45) 





The approximation is refined until the wave number correspond- 
ing to the first mode is determined. The approximation for 
the second mode is twice the value of Ky. The initial esti- 
mate of the remaining eight eigenvalues is the previous 
eigenvalue plus the difference between the previous two 
elgenvalues. 

ievenmineritirst three eigenvalues are determined, the 
program plots their modal structure. The remaining seven 


wave numbers are then listed. 


2 Il 





As the selected frequencies within the inertio-gravity 
band increase in magnitude, the boundary constraints, though 
analytically well-defined, weaken exponentially. The higher 
frequencies are greater than the incremental Vaisala frequency 
for a greater percentage of column depth and it becomes 
increasingly difficult for the numerical program to find 
eigenvalues. RUNK terminates its search for the next eigen- 
values by (1) stopping after 20 refinements by the Newton- 
Raplison iteration process, (2) restricting the magnitude of 


ae em/sec, and (3) stopping if the 


the eigenfunction to 10 
solution converges to a mode other than the one searched for. 
Also, as the frequency within the band increases, W 
becomes more sensitive to changes in k when seeking higher 
moces. For situations where a minute change in k causes 
abrupt responses in end values of W, the eigenvalue may be 


only refined to a certain degree. These eigenvalues are 


indicated as approximations. 


ae 





IV. TESTING THE PROGRAM 


A. LINEAR INPUT 

To check numerical formulation, coding, and accuracy, two 
test cases were used whose analytic solutions were possible. 
The first test was a linear variation of depth from 1.023 
grams/cc at the surface to 1.033 grams/cc at the bottom depth 
of 1000 meters. The Vaisala frequency for this input was 


constant and equal to: 


Ney ae ee 


2 oe 


a Some sec -5 —3 -1 
= 17055 614 mo le =-9,77X10-~sec”. (46) 


The piot by the program of the vems1lty matches that of the 
linear variation. The Vaisala frequency determined numeri- 
cally matches that determined analytically. Table I compares 
the analytical and numerical results of the linear variation 
of density for frequency 0.00684 cycles/min. Figure 4 
compares the first three modal structures. Analytically, the 


normal modes are 


K =( (47) 
m N? — £° h 
tie -ergentunetion 1s: 
2 
NU Z mil 
exp oF i Sala oro Zi) (48) 


B. THREE-LAYER INPUT 
The second input was a density strata that was constant 


at 1.023 grams/cc from the surface to 300 meters, varied 


oe 





Table I. A comparison of the first ten normal modes 


determined for frequency 0.00684 cycles/min: 


Analytical Numerical % rror 
ky 3.613x10> 3.629x 10° 0.45 
k, (BOS sone 7.259x107> 0.47 
k 1.083 x 107" 1.088 x 107 ney 
ky 1.445 x 107 See On 0.42 
ks isco 10  ¢ pil so 0.44 
ke DME 107 2.177 x 107" 0.46 
ke 2.529 x 107 2.540 x 10) 0.43 
kg 2.890x10* A GOS sc 1 0.45 
Kg 3.251 x 107" 3.266 x 107 0.46 
Ki9 Be ols x 10 562 ora Oo 0.44 
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linearly from 300 meters to a value of 1.027 grams/cc at 700 
meters, and then remained invariant from 700 meters to the 
bottom at a depth of 1000 meters. The corresponding value 
of the Vaisala frequency was then zero in the upper and lower 


isopycnal layers and equal to ee = 


in the middle region. 
The density plot by the program matched the analytic 
solution to the third decimal place. Figure 5 illustrates 
the comparison between the analytic solution and the computer 
solution for the Vaisala frequency, and illustrates how the 


interpolation scheme handles discontinuities in the density 


slope. 
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FIGURE 5. Comparison of Vaisala Frequency Plots 
by Analytical and Numerical Methocs 





V. APPLICATION OF THE PROGRAM TO OCEANOGRAPHIC DATA 


The program utilizes density strata from Nansen casts. 
It can process several stations. IDATA, the number of sta- 
tions to be studied, is the first card to be typed. It 
follows the alphanumerics that are read in and Serves in the 
plotting portions of the program. IDATA is read in as an 
integer, right-oriented to column ten of the IBM card. The 
Second card to be typed up consists of three figures. The 
trot 1s M, the number of bottles at the station, read in as 
an integer, right-oriented to column 10. The second is ZMAX, 
the depth at which the last bottle was taken, read in the 
MmomMaver iO. jyerecht-orlente: te column 2G. Vhe vhirde figure 
is RLAT, the latitude of the station, that is read in using 
horiiau FLO, rpreht—oriented to column 30. 

The next cards to be typed are those representing the 
density strata. These are read in as four pairs of numbers 
for each card. Each number has a field width of ten and is 
right-oriented. The first number in each pair represents 
iiiomecdeptieand has a tormat of Fl1O.1. The second of the 
pair has a format F10.5 and represents the station density. 
These four pairs of numbers, representing 80 columns per card, 
are entered until all the station data has been read in. 

The above represent all the cards needed to be typed for 
each station. From experience gained using the IBM 360 


computer to process ocean stations taken in mid-August, 1970, 
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off Monterey, California, it appears that two minutes are 
required to process shallower stations (around 200 meters) 
and four minutes to process deeper stations (around 2000 
meters). 

The actual printed output for each station appears as 
follows. On the first page of printed output appears: 
Gieeethemnumber of station depths, (2) the lowest bottle depth, 
(3) a printout of the station depths and densities in columnar 
form, (4) the inertial frequency computed for the station, 
(5) the maximum Vaisala frequency computed for the station, 
and (6) the ten frequencies within the inertio-gravity band 
selected for consideration. 

The second page of printed output plots the Vaisala 
frequency and density versus dceptn. 

The next ten pages plot the ray path for each frequency 
versus depth. At the bottom of each ray plot, the amount of 
time that the frequency is greater than the Vaisala frequency 
momindica ted. 

The remaining twenty pages of output display information 
on normal modes. Corresponding to each of the ten selected 
frequencies are two pages. The first page lists the first 
three normal modes and plots these modes versus depth. The 
second page lists the remaining seven modes. 

Two stations, one shallow and one deep, were chosen to 
dllustrate the output from the program. fhe results for both 


Slavionseare summarized by the following tables and figures. 
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Table II lists the information for the deep oceanographic 
station displayed on the first page of computer output. This 
station is representative of an oceanic internal wave regime. 
Figure 6 shows the plot of the density variation and the 
Vaisala frequency for the deep station. The Vaisala plot 
follows the pattern of the density profile. 

Pivt@uemiscdepreths the ray paths for the first three select— 
pumuregquencies. The first frequency selected, 0.01371 
cycles/min, never exceeds the incremental Vaisala frequency. 
Thus, theoretically, its solution is always oscillatory and 
it travels horizcntally nearly 30 km. The second frequency 
selected, 0.12969 cycles/min, exceeds the Vaisala frequency 
in the region of the water column below a depth of 9/70 meters. 
Thus it exponentially ‘decays' in this region and its ray 
path over this portion of the water column is vertical. 
Figures 7 and 8 illustrate that, with increasing frequency, 
the horizontal distance covered by the ray path decreases 
and a greater portion is in the exponential region. The ray 
paths for the last four selected frequencies were not shown 
since their ray is nearly vertical from top to bottom. 

fiomem plots the first three modes for frequency 0.013/1 
cycles/min. This was the oniy frequency for which the first 
three mornal modes could be positively identified. Besides 
gaining a better understanding of the actual appearance of 
the modal structures, Figure 9 suggests optimum placement of 
a thermocline follower for mode detection. For example, the 


follower would best detect mode three at a depth of about 
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TABLE II. Data display for deep station. 


Incremental depth Incremental density 


O20 1.02498 
HOG, ie 62505 
PONG. iO 2522 
BOO 1025.65 
SOLAS, 1.02588 
13.20 OZ Gul 

MOON, i 02624 
150.0 1.02646 
200.0 1.02656 
300.0 lp O26 09 
400. 0 1.02681 
DOOn lie G2 162 
OOO. O eO2 al 
800.0 IO 2 30 
1000.0 lO 2 7a 
1200.0 Me O27 50 
1 SOG ee Ie2O27 59 
ZOOO-O leO2 fio 
25 OG Veto 4 
3000.0 kee 7G 
Frequency number Cycles/min 
imcieial 0.00082 

il Oe Ola (a 

2 0.12969 

3 O25 05) 

4 Oe3e 041 

») Ceo le2s 

6 O2645 14 

ff QO. 1 £401 

8 O9O25 7 

9 IA OSL Tes} 
10 i 16059 

iaselimlun Vaisala 1.28946 
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2000 meters. Table III lists the various modes determined 
for the selected frequencies for the deep oceanographic 
station. 

Table IV displays the first page of printed output for 
the shallow station representative of the regime in a large 
bay. Figure 10 illustrates the density variation and the 
Vaisala frequency as a function of depth. As in the deep 
station, the first selected frequency of the shallow station 
has an oscillatory solution throughout most of the water 
column. Its ray path travels horizontally nearly 2500 meters. 
Its exponential region may be seen as a vertical path at a 
distance of 1955 meters. Figure 11 shows the effect of an 
increase in frequency without a considerable increase in the 
Seroecy of exponential decay. ‘The second frequency soicctcca, 
0.15818 cycles/min., has an exponential solution over the 
water column in nearly the same region as the first. However, 
frequency 2 only travels a horizontal distance of 230 meters. 
Figure 12 depicts the decrease in ray path as the frequency 
increases. 

Figure 13 depicts the first three normal modes for 
peequencey O,01656 cycles/min. Figure 14 illustrates the 
first three normal modes for frequency 0.31554 cycles/min. 
Figure 14 illustrates well the manner in which the modal 
structure deteriorates with increased frequency. Table V 
lists the modes determined for selected frequencies at the 


Snewkow Stablion . 
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Table I1lI1. Normal modes for selected frequencies at 


the deep station. 
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Table IV. Data display for shallow station. 


Incremental depth Incremental density 


O20 PO 2535 

Dra e250 f 
LOGO 1.02543 
L><0 1202551 
20.0 O20 

CL Sere) P2029 50 
3050 1.02606 
40.0 iPO26 10 
DOO i O2614 
60.0 1026014 
(35.0 PO 2olko 
100.0 1.02619 
Frequency number Cycles/min 
Tne ia 0.00083 
i C2 Ciecsc 

2 Oley ipsyeul Ke’ 

3 OF 1554 

4 Ora 239 

> 0.63024 

6 OL, eke 

7 0.94495 

8 LOZ 3a 

9 1.25966 

10 ee Oe 
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Table V. Normal modes for selected frequencies at the 


shallow station. 


SeOlooG cycles/min. Ol pcmome ye les/ mini 
ibe 200087 701 1. 0O.0086086 

Qe 0025324 2m 25 oll 
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(Oe 2 36 39 

CeO. eco2oo 
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ieee ©. 31 399 

0.47289 cycles/min, 0.63024 cycles/min. 
1. 0.032306 One SO poo 

peo. 04660 2 Oe Lao 
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Aen 2 O40 Ol s1 60 sev eles min. 
Sees SO L6 

63 CORO veya lO. 0924563 
O45 5.65 2 Oe oe 

oy 4 Dererer4 

9. 0.63999 

Om Ow 2o20 
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VI. CONCLUSIONS 


The numerical program developed in this thesis has derived 

a technique that links the theoretical application of internal 
waves to real data. Real solutions are determined for actual 
density profiles. The program can yield a useful description 
of the internal wave regime within a given area. From 
representative samples of stations investigated, four points 
are apparent: 

(1), that with greater depth the effect of the boundary 

conditions weakens; 

(2). that for low frequencies that exceed the Vaisala 

frequency voveF outy & smail percentage of tno water coluumn, 

the ray path approaches the horizontal, while for higher 

frequencies the ray path approaches the vertical; 

(3), that the plot of the lowest three modes demonstrates 

vieminanimer im which, with increasing frequency, the effect 

of the boundary conditions weakens as the oscillatory 

regions give way to exponential "decay" regions; 

(4), that for the stations investigated by this program, 
only those frequencies within the lower quarter of the 
inertio-gravity frequency band could be solved for normal 


modes with any degree of success. 
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ITS GREATER THAN THE VAISALA FREQUENC 
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